library(tidyverse)
library(ggpubr)
library(DESeq2)

jzl_list <- read_csv("Archive/FPP_TRPV3/data/JZL_list.txt")

fpp_gene <- read_csv("Archive/BD-MassSpec/data/fpp_genes.txt") |>
  pull(1)

care_list <- c(jzl_list[[1]], jzl_list[[2]], fpp_gene)

kls_count <- read_csv("Archive/FPP_TRPV3/data/Thanert-Pieper_2019_NSTI.csv", name_repair = make.names)

kls_meta <- read_delim("Archive/FPP_TRPV3/data/Thanert_meta.txt", name_repair = make.names)

kls_count_df <- kls_count |>
  select(-c(2:3)) |>
  column_to_rownames("Ensemble.gene.ID")

dss_meta <- kls_meta |>
  mutate(
    Sample = make.names(Sample),
    Location = replace_na(Location, "NA") |> make.names()
  ) |>
  select(c(Sample, Location, Classification)) |>
  column_to_rownames("Sample")

dss_meta$Location |> table()
dss_meta$Classification |> table()

multi_dss_meta <- dss_meta |>
  filter(Classification %in% c('Others','Streptococcus'))

multi_kls_count <- kls_count_df |>
  select(rownames(multi_dss_meta))

# force a round integer count matrix -----------
dds <- DESeqDataSetFromMatrix(
  countData = round(multi_kls_count),
  colData = multi_dss_meta,
  design = ~ Location + Classification
)

# pre-filtering data of too low counts
keep <- rowSums(counts(dds)) > 1
dds <- dds[keep, ]

# relevel to specify the control group (reference)
levels(dds$Classification)
dds$Classification <- relevel(dds$Classification, ref = "Others")

# rlog() may take a few minutes with 30 or more samples,
# vst() is a much faster transformation
vsd <- vst(object = dds, blind = TRUE)

dds <- DESeq(dds)

res <- results(dds)
summary(res)

normalized_counts <- assay(vsd)

as.data.frame(res) %>%
  rownames_to_column("Ensemble.gene.ID") %>%
  left_join(kls_count[1:2]) -> deg_tibble

deg_tibble |>
  filter(Official.gene.symbol %in% care_list & pvalue < 0.05)

deg_tibble |>
  ggplot(aes(log2FoldChange, -log10(pvalue))) +
  geom_point()

write_csv(deg_tibble, "Archive/FPP_TRPV3/results/")

deg_tibble %>%
  mutate(type = case_when(
    log2FoldChange > 2 & padj < 0.05 ~ "Enriched in infection",
    log2FoldChange < -2 & padj < 0.05 ~ "Depleted in infection",
    .default = "NS"
  )) -> deg_tibble

ggplot(deg_tibble, aes(log2FoldChange, -log10(pvalue))) +
  geom_point(aes(color = type)) +
  geom_vline(xintercept = c(2, -2), linetype = "dashed") +
  geom_hline(yintercept = 3.5, linetype = "dashed") +
  scale_color_manual(values = c("red", "blue", "grey"))

# clusterProfiler -----------
deg_tibble %>%
  filter(log2FoldChange > 0, padj < 0.1) %>%
  pull(SYMBOL) -> up_list

deg_tibble %>%
  filter(log2FoldChange < 0, padj < 0.1) %>%
  pull(SYMBOL) -> down_list

enrichGO(down_list,
  OrgDb = "org.Hs.eg.db",
  keyType = "SYMBOL",
  ont = "ALL",
  universe = deg_tibble$SYMBOL,
  readable = TRUE
) -> ora_go

dotplot(ora_go)

cnetplot(ora_go)

enrichplot::pairwise_termsim(ora_go) %>%
  enrichplot::treeplot()
